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Abstract 

Currently, the high-precision estimation of nonlinear parameters such as Gini in- 
dices, low-income proportions or other measures of inequality is particularly crucial. 
In the present paper, we propose a general class of estimators for such parameters that 
take into account univariate auxiliary information assumed to be known for every unit 
in the population. Through a nonparametric model-assisted approach, we construct a 
unique system of survey weights that can be used to estimate any nonlinear parameter 
associated with any study variable of the survey, using a plug-in principle. Based on a 
rigorous functional approach and a linearization principle, the asymptotic variance of 
the proposed estimators is derived, and variance estimators are shown to be consistent 
under mild assumptions. The theory is fully detailed for penalized B-spline estimators 
together with suggestions for practical implementation and guidelines for choosing the 
smoothing parameters. The validity of the method is demonstrated on data extracted 
from the French Labor Force Survey. Point and confidence intervals estimation for 
the Gini index and the low-income proportion are derived. Theoretical and empirical 
results highlight our interest in using a nonparametric approach versus a parametric 
one when estimating nonlinear parameters in the presence of auxiliary information. 
Keywords auxiliary information; penalized B-splines; calibration; concentration and 
inequality measures; influence function; linearization; model-assisted approach; total 
variation distance. 
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1 Introduction 



The estimation of nonlinear parameters in finite populations has become a crucial problem 
in many recent surveys. For example, in the European Statistics on Income and Living 
Conditions (EU-S1LC) survey, several indicators for studying social inequalities and poverty 
are considered; these include the Gini index, the at-risk-of-poverty rate, the quintile share 
ratio and the low-income proportion. Thus, deriving estimators and confidence intervals 
for such indicators is particularly useful. In the present paper, assuming that we have a 
single continuous auxiliary variable available for every unit in the population, we propose 
a general class of estimators that take into account the auxiliary variable, and we derive 
their asymptotic properties for general survey designs. The class of estimators we propose 
is based on a nonparametric model-assisted approach. Interestingly, the estimators can be 
written as a weighted sum of the sampled observations, allowing a unique weight variable 
that can be used to estimate any complex parameter associated with any study variable of 
the survey. Having a unique system of weights is very important in multipurpose surveys 
such as the EU-SILC survey. 

The estimation of nonlinear parameters is a problem that has already been addressed 
in several papers such as Shao (1994) for L-estimators, Binder and Kovacevic (1995) for 
the Gini index and Berger and Skinner (2003) for the low-income proportion. We mention 
also the very recent work of Opsomer and Wang (2011). Taking auxiliary information into 
account for estimating means or totals is a topic that has been extensively studied in the 
literature; it now encompasses the model-assisted and the calibration approaches, which 
coincide in particular cases (Sarndal, 2007). In a model-assisted setting, linear models are 
usually used, thus leading to the well-known generalized regression estimators (GREG). Some 
nonparametric models have also been considered (Breidt and Opsomer, 2009). However, to 
the best of our knowledge, ratios, distribution functions and quantiles are the only examples 
of nonlinear parameters estimated using auxiliary information. 

To derive our class of estimators and their asymptotic properties, we use an approach 
based on the influence function developed by Deville (1999). This approach utilizes a func- 
tional interpretation of the parameter of interest and a linearization principle to derive 
asymptotic approximations of the estimators. In general, the precision of an estimator $ 
of a nonlinear finite population parameter $ is obtained by resampling techniques or lin- 
earization approaches and in the present paper we focus on linearization techniques. When 
a sample s is selected from the finite population U according to a sampling design p(-), the 
linearization of $ leads under some assumptions, to the following approximation: 



where = Pr(k e s) > denotes the first-order inclusion probability for element k under 
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the design p(-). The right term of ([T]) is the difference between the well-known Horvitz- 
Thompson estimator and the parameter it estimates, namely the total of the variable Uk 
over the population U. Here, Uk referred to as the linearized variable of $ and the way 
it is derived depends on the type of linearization method used which could include the 
Taylor series (Sarndal et al, 1992), estimating equations (Binder, 1983) or influence function 
(Deville, 1999) approaches. The artificial variable Uk is used to compute the approximative 
variance of $ as 

VV(tth-™) — — , (2) 

^— ' ^ 7T fe Til 

s s 

with ttm = Pr(k 6 s, I £ s) the joint inclusion probability for the elements k,l £ U. 

Roughly speaking, when examining Q and @, we can see that, if we estimate in an 
efficient way Ylu u k-> we wm acn i eve a small approximative variance and good precision for 
$. As stated above, it is well known that auxiliary information is useful for improving on 
the estimation of a total in terms of efficiency and, based on a linear model, the use of a 
GREG estimator is the most common alternative. When estimating a total, note that the 
asymptotic variance of the GREG estimator depends on the residuals of the study variable 
on the auxiliary variable. Because linearized variables may have complicated mathematical 
expressions, fitting a linear model onto linearized variables may not be the most appropriate 
choice. This may occur even if the study and the auxiliary variables have a clear linear rela- 
tionship, as illustrated in the following example. Consider a data set of size 1000 extracted 
from the French Labor Force Survey and consider (the wages of person k in 2000) as the 
study variable and x^ (the wages of person k in 1999) as the auxiliary variable. We now 
consider the problem of estimating the Gini index. The expression of the linearized variable 
Uk, k £ U for the Gini index is given in Binder and Kovacevic (1995) and recalled in equa- 
tion (17). It is a complex function of the study variable y^, k & U . In the left (resp. right) 
graphic of Figure [TJ the study variable y^ is plotted (resp. the linearized variable Uk) on the 
y-axis and the auxiliary variable Xk is plotted on the x-axis. The relationship between the 
study variable and the auxiliary variable is almost linear; however the relationship between 
the linearized variable of the Gini index and the auxiliary information is no longer linear. 
The consequence of this is that we cannot increase the efficiency of estimating a Gini index 
if we take the auxiliary information into account through a GREG estimator. Therefore, 
nonparametric models should be preferred to estimate nonlinear parameters Recent work 
already employs nonparametric models to estimate totals (Breidt and Opsomer, 2000, Breidt 
et al., 2005 and Goga, 2005). The use of nonparametrics prevents model failure; however 
the improvement over parametric estimation for totals and means may not be significant 
enough to justify the supplemental difficulties of implementing nonparametric methodol- 
ogy. As illustrated above, the motivation for using nonparametrics becomes much stronger 
when estimating nonlinear parameters. Note that the use of nonparametric regression to 
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Figure 1: Left plot: the wages of person k in 2000 against Zk- the wages of person k in 
1999. Right plot: Uk- linearized variable of the Gini index for the wages in 2000 for person 
k against Zk- the wages of person k in 1999. 

estimate distribution functions and quantiles has also been studied, for example in Johnson 
et al. (2008); however, to our knowledge, this has not been performed for other nonlinear 
parameters. 

We propose a novel methodology that allows for the efficient estimation of any parameter 
<3> by combining the functional approach (Deville, 1999) with any of the previously suggested 
nonparametric methods. One issue with the functional approach is that several technical 
details are not provided in Deville (1999); thus it is difficult to derive rigorous proof of asymp- 
totic results by following this approach. In the present paper, we propose to clarify some 
important points and derive rigorous proofs of our asymptotic results. Most importantly, we 
prove that the total variation distance between finite measures is an adequate choice for the 
derivation of asymptotic approximations in this context. Asymptotic results are detailed at 
length for penalized B-spline nonparametric estimators. 

The estimators under study combine two types of nonlinearity: nonlinearity due to the 
expression of a complex parameter and nonlinearity due to nonparametric estimation. We 
propose a two-step linearization procedure that provides an approximation of the nonpara- 
metric estimator via a Horvitz-Thompson estimator of a total using an artificial variable. 
Roughly speaking, this artificial variable corresponds to the residuals of the linearized vari- 
able Uk on the fitted values under the model. Because the linearized variables depend on the 
parameter of interest, the residuals will also depend on this parameter. The consequence of 
this important and general property is that the nonparametric approach helps to get a unique 
system of weights that may lead to a gain in efficiency for different complex parameters. 

The paper is structured as follows: the second section provides some background informa- 
tion on the nonparametric estimation of a finite population total in a general framework. In 



the third section, a class of nonparametric substitution estimators based on nonparametric 
regression is introduced. Variance approximations are derived using the influence function 
linearization approach (Deville, 1999) in a general nonparametric setting. We propose in 
the fourth section a penalized B-spline model-assisted estimator for the finite population 
totals which is in fact an extension to a survey sampling framework of the penalized B-spline 
estimator studied in Claeskens et al. (2009). We prove that the estimator is asymptot- 
ically design-unbiased and consistent. Next, we build the nonparametric penalized spline 
estimation for nonlinear parameters and we assess the validity of the two-step linearization 
technique. The fifth section defines a class of consistent variance estimators while section 
six contains a case study. The data set is extracted from the French Labor Force surveys 
of 1999 and 2000 as presented previously. Asymptotic and finite-sample properties of the 
regression B-spline estimators are illustrated for the simple random sampling without re- 
placement and the stratified simple random sampling. This section also includes suggestions 
for practical implementation and guidelines for choosing the smoothing parameters. Finally, 
section seven concludes this study and the assumptions and the technical proofs together 
with some discussion are provided in the Appendix. 

2 Nonparametric model-assisted estimation of finite 
population totals 

We focus on the estimation of the total t y = Ylk=i Vk = Tlu V k °f the study variable y over 
U, taking into account the univariate auxiliary variable Z. The values zi, . . . , zn of Z are 
assumed to be known for the entire population. 

Many approaches can be used to take into account auxiliary information Z and thus im- 
prove on the Horvitz-Thompson estimator t Vj HT — Yls Uk/^k- The goal is to derive a weighted 
linear estimator t yw = J2 s w ksyk of t y , such that the sample weights w ks do not depend on 
the study variable values y k but include the values z k) for all k G U. The construction of the 
model-assisted (MA) class of estimators i yw is based on a superpopulation model £: 

f: Vk = f\z k ) + e k (3) 

where the e k are independent random variables with mean zero and variance v(z k ). If f(z k ) 
was known for all k G U, the total t y may be estimated by the generalized difference estimator 
(Cassel et al, 1976), 

Note that i yjdiS consists in the difference between the Horvitz-Thompson estimator t Vj HT an d 
its bias under the model £, namely Y^ s f(z k )/n k — f(z k ). As a consequence, t yjdiS is 
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unbiased under the model, E^(t y>dm ) = t y and moreover, it is unbiased under the sampling 
design, E p (t y;diB ) = t y . The variance of i yAm under the sampling design is given by 

t//+ ^ V" \^ ^ Vk - f(z k )yi - f(zi) , * 

^ — ' z — ' 7TA- 7T; 

[/ [7 ret 

which shows clearly that the difference estimator t y>dig is more efficient than the Horvitz- 
Thompson estimator t y> HT if f(zk) approximates well yt for all k & U. 

In practice, we don't know the true regression function /, thus we use an estimator of it. 
Generally, this estimator is obtained using a two-step procedure: we estimate first / by / 
under the model £ and next, we estimate / by / using the sampling design. Plugging / in 
ph, yields the final estimator of t y . 

The linear regression function f(z] C ) = z' k /3 yields the generalized regression estimator 
(GREG) extensively studied by Sarndal et al. (1992). The GREG estimator is efficient if 
the model fits the data well, but if the model is misspecified, the GREG estimator exhibits no 
improvement over the Horvitz-Thompson estimator and may even lead to a loss of efficiency. 
One way of guarding against model failure is to use nonparametric regression which does 
not require a predefined parametric mathematical expression for /. 

Recently, Breidt and Opsomer (2000) proposed local linear estimators and Breidt et al. 
(2005) and Goga (2005) used nonparametric spline regression. The unknown / function is 
approximated by the projection of the population vector = (yi, . . . ,7/at)' onto different 
basis functions, such as the basis of truncated gth degree polynomials in Breidt et al. (2005) 
and the B-spline basis in Goga (2005). In the following, we briefly recall the definition 
and the main asymptotic properties of nonparametric model-assisted estimators for finite 
population totals (see also Breidt and Opsomer, 2009). 

Let f Vi k be the estimator of f(zk) obtained at the population level using one of the three 
nonparametric methods mentioned above. Plugging f y ^ into Q results in the following 
nonparametric generalized difference pseudo-estimator of the finite population total: 

<U = £^^ + £/»> w 

k u 

Note that t* diff is called a pseudo-estimator because it is not feasible in practice since f y ^ 
is unknown. This pseudo-estimator is still design-unbiased but it is model-biased because 
nonparametric estimators f y ^ are biased for f(zk) (Sarda and Vieu, 2000). Nevertheless, 
under supplementary assumptions (Breidt and Opsomer, 2000 and Goga, 2005), the bias 
under the model vanishes asymptotically to zero when the population and the sample sizes 
go to infinity. The unknown quantities f y> k are usually obtained by least squares methods 
(ordinary, weighted or penalized) and we may write 

f y ,k = t&yu, for all k 6 U (7) 
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where the N dimensional vector depends on the population values Zk, k G U as well as 
on the projection matrix for the considered basis functions, but does not depend on y. The 
expression of depends on the chosen nonparametric method, as discussed in Breidt and 
Opsomer (2000), Breidt et al. (2005) and Goga (2005). 

As in the parametric case, we estimate f y ^ by f y ^ using the sampling design, 

fy,k = qtksYs, for all k E U (8) 

where q' ks is the n-dimensional design-based estimator of q' k and y s = (yk)k^s is the sample 
restriction of yu- Plugging f y ^ into (JgJ) yields the following nonparametric model-assisted 
estimator (NMA) 

ty,np = ~~ ~ l~ fy,k- (9) 

s k U 

This estimator can be written as a weighted sum of the sampled observations 

iy,n P = ^2 WksVk = w « ys ' ( 10 ) 

s 

where the weights w s = (wk s )kes depend only on the sample and on the auxiliary information, 



w, 



n^i.-Q'n-^. + Q' iu, (n; 



with l s the n dimensional vector of ones, II S the nxn diagonal matrix with n^, k e s, along 
the diagonal and Qu the Nxn matrix having q' fcs as rows with sample restriction Q s = 



(q' A;s )fc es . The estimator (10) is a nonlinear function of Horvitz-Thompson estimators, and 
its asymptotic variance has been obtained on a case-by-case study. Under mild hypothesis 
(Breidt and Opsomer, 2000, Breidt et al, 2005 and Goga, 2005), i y ^ np is asymptotically 
design-unbiased, namely lim A r_ 5 . 00 -E p (tj /iri p — t y )/N = and design i/n-consistent in the sense 
that 

N-\i y , np -t y ) = O^n- 1 ' 2 ). (12) 

Moreover, it can be approximated by the nonparametric generalized difference estimator 
t* 

N-\i y<np -t y ) = iV-i^-g + o^n- 1 / 2 ). (13) 

Furthermore, if the asymptotic distribution of V^(t* jdiff ) _1 / 2 (t* idiff — t y ) is normal A/"(0,1), 
we have that the asymptotic distribution of V p (t* ytdia )~ l l 2 {t y ^ np — t y ) is also normal ftf(Q, 1) 
where V p (t* Aifi ) is obtained according to formula |5) applied to residuals — f y ^- This 
means that the NMA estimators bring an improvement over parametric methods and the 
Horvitz-Thompson estimator when the relation between y and Z is not linear. In this case, 
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the residuals — f Vt k will be smaller than under a parametric smoother, which explains 
the diminution of the design variance of NMA estimators. Nevertheless, nonparametric 
estimators require that the auxiliary information should be known on the whole population 
unlike the GREG estimator that requires only the finite population total for Z. 

The efficiency of NMA estimators depends on the choice of the smoothing parameters. 
Opsomer and Miller (2005) and Harms and Duchesne (2010) derive the optimal bandwidth 
for the local polynomial regression, while Breidt et al. (2005) circumvent the issue of the 
number of knots by introducing a penalty coefficient. They also give a practical method for 
estimating this penalty. 

3 Nonparametric model-assisted estimation 
of nonlinear finite population parameters 

3.1 Definition of the nonparametric substitution estimator 

Let us consider the estimation of some nonlinear parameters $ by taking into account uni- 
variate auxiliary information known for all the population units. Examples of a nonlinear 
parameter of interest $ include the ratio, the Gini coefficient and the low-income propor- 
tion. A parameter $ may depend on one or several variables of interest; however, the same 
auxiliary variable Z will be used to explain these variables of interest. 
We aim to provide a general method for the estimation of $ using Z and considering the 
functional approach introduced by Deville (1999). The methodology consists in considering 
a discrete and finite measure M = J2u^Vk where 5 Vk is the Dirac measure at the point 
and M is such that there is unity mass on each point y k with k E U and zero mass elsewhere. 
Furthermore, we write $ as a functional T of M, 

$ = T(M). (14) 

The nonparametric weights Wk s are provided by ( [IT] ) and M is estimated by 

M np = ^2 W ks5y k ■ 
S 

Even if these weights are derived to estimate the total t y , they do not depend on the study 
variable y-, thus they can be used to estimate any nonlinear parameter of interest $ when 
it can be expressed as a function of M. Note that M np is a random measure of total mass 
equal to N np = ^ s w ks . 

Plugging M np into (14) provides the following nonparametric substitution estimator for $, 



= T(M; 



np ) 
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We will now illustrate the computation of <3> np using the simple case of a ratio R and subse- 
quently the more intricate case of the Gini index and parameters defined by implicit equa- 
tions. 



a. The ratio R between two finite population totals. We write R = J2u V k l Xk m a 

( ydMiy) 

functional form as R = j — The nonparametric estimator of R is easily obtained 

by replacing the measure M with M np , namely R np = ^ ^ \ np ^ = ^ s Wfcs ^ fc j± s j m ii ar 

/ xdM np (x) 2^ s w ks x k 

estimation of R using GREG weights was previously considered by Sarndal et al. (1992). 

b. The Gini index. The Gini index (Nygard and Sandstrom, 1985) is given by 

c _ EuVk (2F(y k ) - 1) _ J(2F(y) - l)ydM(y) 
t y JydM(y) 

where F(y) = J l{£< y ydM(£) / J dM(y) = J2u ^-{yk<y}/N is the empirical distribution func- 
tion. Again, the nonparametric estimator for G is obtained by simply replacing M with 
M np . Hence, 

X J2s w ks{2F np (y k ) - l)y k 

^np — , [i-O) 

L s w k S yk 

where F n Jy) = [J^^M = Lg! 1 !^, 
JdM np (y) l^ s w ks 

c. Parameters defined by an implicit equation. Let $ be defined as the unique solution 
of an implicit estimating equation X^i/0fe( < ^ ) ) = (Binder, 1983) that may be written 
in a functional form as j 0($)<iM = 0. We replace M with M np and the nonparametric 
sample-based estimator of $ is the unique solution of the sample-based estimating equation 
J <fi{§)dM np = Y2 s Wks< l )k (^ n p) = 0- ^ n exam ple of such a parameter is the odds-ratio which 
is extensively used in epidemiological studies. Goga and Ruiz-Gazen (2012) have studied the 
estimation of the odds-ratio by taking into account auxiliary information and nonparametric 
regression. 



3.2 Asymptotic properties of the nonparametric substitution es- 
timator under the sampling design 

In this section, we investigate the asymptotic properties of the nonparametric estimator 
using the asymptotic framework suggested by Isaki and Fuller (1982). Additionally, we make 
several assumptions (detailed in the Appendix) regarding the regularity of the functional T 
and the first order inclusion probabilities of the sampling design. 
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The nonparametric estimator $ np is doubly nonlinear, with nonlinearity due to the pa- 
rameter $ and nonlinearity due to the nonparametric estimation. Our main goal is to 
approximate $ np using a linear estimator (Horvitz-Thompson type) which will allow to com- 
pute the asymptotic variance of This approximation will be accomplished in two steps: 
first, we will linearize $ and next, we will linearize the nonparametric estimator obtained in 
step one. 

The first linearization step is a first-order expansion of $ np with the reminder going to 
zero. The parameter of interest $ is a statistical functional T defined with respect to the 
measure M or equivalently, with respect to the probability measure M/N (by assumption 
Al). Using the first-order expansion of statistical functionals T as introduced by von Mises 
(1947) and under the assumption of Frechet differentiability of T, the reminder depends on 
some distance function between M/N and an estimator of this measure (Huber, 1981). Dev- 
ille (1999) uses these facts to prove the linearization of the Horvitz-Thompson substitution 
estimator of $; however, no details are given about the considered distance, while Goga et 
al. (2009) provide only minimal details. In what follows, we provide a distance between 
M np /N and the true M/N which goes to zero when the sample and the population sizes go 
to infinity. 

We consider the total variation distance for two finite and positive measures Mi and M2 
to be defined by 



distance d^ v between the Horvitz-Thompson estimator of M/N and the true M/N goes to 
zero. Next, we extend the result (lemma [2] from below) to the nonparametric estimator 



Let Wks represent the Horvitz-Thompson weights, namely Wk s = 1/iik for all k G s and let 
Mjjx = J2 s ^yk/ n k De the estimator of M using these weights. Let h G H and for ease of 
notation, Xk = h(yk). Thus, for all k G U, \xk\ < 1 uniformly in h G % and 




with % = {h : R — > R| sup^ \h(x)\ < 1}. We first prove (lemma 



I] from below), that the 



M np /N. 



I 



hdM, 



HT — 



I 




where Ik = l{fc £s } is the sample membership indicator. 



Lemma 1. Assume (A3) and (A5) from the Appendix. Then, 




The proof is provided in the Appendix. We extend now lemma [T] to nonparametric weights 
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Wk s given by (11). Consider again h eH and let 



J h dM np - J hdM = 



u 



+ 2^/^-2^ Xk 

s k U U 

where is obtained from (Ji| for y k replaced with x k = h{y k ). Let also f Xjk obtained from 
([7]) for y k replaced with x k . 

Lemma 2. Assume (A3) and (A5) from the Appendix. Assume in addition that: 
(A*) for all k e U, ^ J2u fx,k = uniformly in h. 

(A**) E p |i Ec/ (| - X ) (/«,* " = Ofa- 1/2 ) imi/orrofy zn fc. 
T/ien, 

d to (m^p/JV, M/iv) = O^n- 1 / 2 ). 

The proof is provided in the Appendix. In section |4j we prove that the nonparametric 
estimator of M constructed using B-spline estimators satisfies the assumptions (A*) and 
(A**) from the above lemma. The results from Breidt and Opsomer (2000) may be used 
to prove the assumptions for local polynomial regression; however, this issue will not be 
pursued further here. 

To provide the first order expansion of $ = T(M), we must also define its first derivative. 
This derivative is referred to as the influence function and is defined as follows (Deville, 1999) 

where 5 y is the Dirac measure at point y. Note that the above definition is slightly different 
from the definition of the influence function given by Hampel (1974) in robust statistics, 
which is based on a probability distribution instead of a finite measure. 
Let Uk, for all k £ U be the influence function IT computed at y — y k , namely 

u k = IT{M,y k ), keU. 

These quantities are referred to as the linearized variables of $ and serve as a tool for 
computing the approximative variance of $ np . They depend on the parameter of interest 
and they are usually unknown even for the sampled individuals. Deville (1999) provides 
many practical rules for computing u k for rather complicated parameters 
Examples. The linearized variable of a ratio R is 

u k = =^ — (y k -Rx k ) (16) 
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and for the Gini index, it is given by 

, y k -y k ,< l + g , i-G 

u k = 2F(y k ) y k — — + (17) 

where y k < is the mean of yj lower than y k . 

We now provide the main result of this paper. The following theorem is the first lineariza- 
tion step of This proves that under broad assumptions the nonparametric estimator 
$ np is approximated by the nonparametric estimator for the population total J2u u k °f the 
linearized variable. The proof is provided in the Appendix. 

Theorem 3. (First linearization step) Assume (A1)-(A3) and (A5) from the Appendix. 
Additionally assume (A*) and (A**) from lemma^ Then, the nonparametric substitution 
estimator $ np fulfills 

N~ a ($ np -^ = N~ a n^w ks u k -J2 u k) + o P (n~ 1/2 ). 
We can put J2 S w ks u k in the form of an NMA estimator. Let denote £* = w ks u k . Using 



(11), we can write 



u,np - w s n s — / , r } y g u £, 

s k U 



t.. „„ = w„ 



where g* u k = q' fes u s with q ks is given by (JSJ) and u s = (u k ) k&s is the sample restriction of 

Remark 1: A model-based interpretation of g* uk may be given. For the nonparametric model 
the linearized variable u k can be fitted using the auxiliary variable z k , 

f ' : u k = g{z k ) + 7] k 

where the r] k are independent random variables with mean zero and variance v(z k ). The 
estimator of g under the model denoted by g U)k , is obtained using the same nonparametric 
method employed for estimating f under the model £. This implies that g U)k = c\' k uu is the 
best fit of the population vector U[/ = (u k )keu with q k given by ([?]). Furthermore, q k is 
estimated by q ks which leads to the pseudo-estimator = q' fcs u s of g Ujk . However, unlike 
the linear case, g* k is not an estimate of g Uyk because the sample linearized variable vector 
u s is not known and we refer to it as a pseudo-estimator. Remark also that the estimator 
$ np is efficient if the nonparametric model £' holds. 



The nonparametric pseudo-estimator t* unp given by (18) is a nonlinear function of Horvitz- 
Thompson estimators; however, it estimates a linear parameter of interest, namely the total 
of u k , t u = J^jjUk. This indicates that t* unp is similar to estimators used by Breidt and 
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Opsomer (2000), Breidt et al. (2005) and Goga (2005) although it is computed for the 
artificial variable Uk- The second linearization step approximates t* unp by the generalized 
difference estimator of Ylu Uk gi ven by 

t:, diS = + ]>>,,. ( 19 ) 

s k U 

Proposition 4. (Second linearization step) Assume that N~ a {t* u ^ np — t^ diff ) = o^n' 1 / 2 ). 
Then, 

N- a (C np - t u ) = N- a (t* - t u ) + o^n- 1 ' 2 ). 



Based on theorem |3j and proposition 4 , we see that the asymptotic variance of $ np is the 
variance of t* u d;ff , namely 



u k - g u ,k ui - g U)i 

kV 



7Tk 7T; 

U U 

Moreover, if the asymptotic distribution of Vp 2 (t^ diff )(£„ diff — t u ) is jV(0, 1), then the asymp- 

totic distribution of V p (*ud»)(^p — ^) * s a ^ so -^"(0) !)■ m section 4, we provide the nec- 
essary assumptions for the linearized variables and the auxiliary variable Z to obtain an 
approximation of t* u np by £jj jdiff in a B-spline estimation context. 

Remark 2. When the linearized variable Uk is a linear combination of the study variables, 
the assumption from proposition [7] is reduced to assumptions on the study variables. For 
example, this occurs in the case of a ratio R = t y /t x , where the linearized variable is given 

by Uk = —(Vk — Rxk) = Aiy k + A 2 x k . The error t^ diff — £* n can be written as a linear com- 

bination of errors between t* idiff — ty^ np and t x>diff — t x>np , respectively. Using mild regularity 
assumptions on X, y and on the sampling design, N^ity^p — t* diff ) and N~ l {t x ^ np — t xdiff ) 
are shown to be of order o p (n -1 / 2 ) (see Fuller, 2009, for linear regression and section^ for 
B-spline estimators) . Thus t* u — t* u diff is also of order o p (n~ l l 2 ) provided that R and N~ 1 t x 
are bounded. 

Remark 3. The asymptotic variance Q np given by theorem^ and proposition 4 depends on 
the population residuals Uk — g u ,k of the linearized variables Uk under the model For the 
simple case of a ratio, the relationship between Uk and the study variables is explicit and given 
by Uk = Aiyk + A 2 Xk. If linear models fit the data Xk and yk well, then a linear model will also 
fit Uk well. Nevertheless, for nonlinear parameters such as the Gini index, the relationship 
between Uk and the study variable is not as simple as that for the ratio. In such situations, 
the use of nonparametric regression methods may provide a major improvement with respect 
to variance compared to parametric regression. 
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4 Penalized B-spline estimators 



Spline functions have many attractive properties, and they are often used in practice due 
to their good numerical features and ease of implementation. We suppose without loss of 
generality that all z^ have been normalized and lie in [0, 1]. For a fixed m > 1, the set Sk,™, 
of spline functions of order m, with K equidistant interiors knots = £o < £i < • • • < Ck < 
= 1 is the set of piecewise polynomials of degree m — 1 that are smoothly connected 
at the knots (Zhou et ai, 1998), 

Sr,™ = {t € C m ~ 2 [0, 1] : t(z) is a polynomial of degree (m— 1) on each interval [£i,£t+i]} 

For m = l, Sx,m is the set of step functions with jumps at knots. For each fixed set of knots, 
Sk,™ is a linear space of functions of dimension q = K + m. A basis for this linear space is 
provided by the B-spline functions (Schumaker, 1981, Dierckx, 1993) Bi, . . . , B q defined by 

t=0 11 r=0,r^l\ c ij-l Sj-r) 

where — x)™ ~ 1 = — x) m_1 if > x and zero, otherwise. For all j = 1, . . . , q, 
each function Bj has the knots ^_ m , . . . , ^ with ^ r = £ min ( m£LX (r,o),K+i) for r = j - m, . . . ,j 
(Zhou et ai, 1998) which means that its support consists of a small, fixed, finite number of 
intervals between knots. Moreover, B-spline are positive functions with a total sum equal to 
unity: 

J2b j (x) = i, xe[o,i}. (20) 

i=i 

For the same order m and the same knot location, one can use the truncated power basis 
(Ruppert and Carroll, 2000) given by 1, z,z 2 ,..., z m ~\ (z - £1)+^, ...,{z- + The 
B-spline and the truncated power bases are equivalent in the sense that they span the same 
set of spline functions Sk™ (Dierckx, 1993). Nevertheless, as indicated by Rupert et al. 
(2003), "the truncated power bases have the practical disadvantage that they are far from 
orthogonal" , which leads to numerical instability especially if a large number of knots are 
used. 

4.1 Nonparametric penalized spline estimation for finite popula- 
tion totals 

We now consider the superpopulation model £ given by (|3j). To estimate the regression 
function /, we use spline approximation and a penalized least squares criterion. We define 
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the spline basis vector of dimension q x 1 as b'(zk) = (Bi(z k ), ■ ■ ■ , B q (z k )), k G U. The 
penalized spline estimator f y>k of f(z k ) is given by f y ^ k = b'(z k )Q y) \ with yX as the least 
squares minimizer of 

N x 

yi Vk - h'{z k )ef + A / [(W(t)6)^} 2 dt, (21) 
fc =i Jo 



where ^ represents the p-th derivate with p < m — 1. The solution of (21) is a ridge-type 
estimator, 

y>A = (B'^Bf/ + XD^B'uyu, (22) 

where B[/ is the N x q matrix with rows b'(z k ) and the q x q matrix D p is the squared 
L 2 norm applied to the pth derivative of h'O. Because the derivative of a 5-spline function 
of order m may be written as a linear combination of S-spline functions of order m — 1, 
for equidistant knots D p = K 2p VpRV p (Claeskens et ai, 2009) where the matrix R has 
elements Rij = B^" 1 '^ (t) Bj m ~ p \t)dt with B^ m ~ p ^ as the i?-spline function of order m — p 
and V p as the matrix corresponding to the pth order difference operator. 

The amount of smoothing is controlled by A > 0. The case A = results in an unpenalized 
B-spline estimator the asymptotic properties of which have been extensively studied in the 
literature (Agarwal and Studden, 1980, Burman, 1991, and Zhou et ai, 1998, among others). 
The case A — » oo is equivalent to fitting a (p — l)th degree polynomial. The theoretical 
properties of penalized splines with A > 0, have been studied only recently by Cardot (2000), 
Hall and Opsomer (2005), Kauermann et al. (2009) and Claeskens et al. (2009). 
The design-based estimators of f y ^ are 

f Vt k = b'(z k )6 y ,x (23) 

where y ^ = (B / S II~ 1 B S + AD p ) _1 B / s Il7 1 ys is the design-based estimator of y X and B s is 
the n x q matrix given by B s = (b'(Zk))kes- We note that f Vik may be written as in formula 
g for q' fcs = b'(^)(B' s nj 1 B s + ad^-^iij 1 . 
Finally, the -B-spline NMA estimator of t y is as follows: 



1 v,bs 



(24) 

This indicates that t Vt BS may be written as a GREG estimator that uses the vectors b'(z k ) 
as regressors of dimension q x 1 with q going to infinity and a ridge-type regression coefficient 
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O y ,\- Furthermore, t Vt BS is a weighted sum of sampled values y k with weights w s expressed 
as in (11), 



w s = n; 1 !, - ii^b, (Bln^B, + ad p ) 1 (B'.n; 1 !, - b'^) 



(25) 



Regression splines 

For A = 0, we obtained the unpenalized B-spline estimator studied by Goga (2005) and called 



the regression splines. The B-spline property given in (20) may be written as V ■ h(z k ) = 1 
with l q the q dimensional vector of ones, implying that l s = B s l g and lu = Bj/lg. Using 



these two relations in (25 ) (Goga, 2005), we observe that t Vt BS is equal to the finite population 
total of the prediction f y>k = b / (^)(B;n; 1 B s )- 1 B / s n7 1 y s 



's, 



iy,BS - ^2 fv> k ~ W '^ s 
U 

where the weights are given by, 

Ws = n; 1 B s (Bln^B,) -1 b^. (26) 

Note the similarity with the GREG weights obtained in the case of a linear model when 
the variance of errors is linearly related to the auxiliary information (Sarndal, 1980). We note 
that for a B-spline of order m = 1, the estimator t Vt BS becomes the well-known poststratified 
estimator (Sarndal et ai, 1992). 

Based on assumptions regarding the sampling design and the variable y, (assumptions (A3)- 
(A5) from the Appendix) and assumptions regarding the distribution of Z and the knot 
number (assumptions (B1)-(B2) in the Appendix), Goga (2005) proved that the B-spline 



estimator for the total t y is asymptotically design-unbiased and consistent (equation (12)) 



and may be approximated by a nonparametric generalized difference estimator (equation 



(13)). These results are valid without supplementary assumptions regarding the smoothness 



of the regression function /. 



Penalized splines using truncated polynomial basis functions 

Let c'(z k ) = {1, z k , . . . , z™' 1 , (zfc-£i)+~\ • • • , 0zfc-6s-)+ _1 } be the vector basis and let f y>k = 
c'(z k )fjy,p with fj y>p be the least squares minimizer of J2k=i(Vk ~ c '( z k)vY + PJ2f=iVm-i+j 
for T]' = (r] , . . . , rj m _i +K ). The solution is given by 

with Cjj = (c'(zk))keu an d the penalty matrix A having m — 1 zeros on the diagonal fol- 
lowed by K one values, A = diag(0, . . . , 0, 1, . . . , 1). Note that for p = 0, we obtain the 
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same prediction f y ^ as with an unpenalized B-spline estimation. This results follows from 
the fact that the two bases are equivalent, thus there exists a square and invertible tran- 
sition matrix Ljj such that B[/ = C{/L[/ (Ruppert et al, 2003). For p > 0, we have 
f Vt k = B' [/ (B' [/ B[/ + pL' u AL u )~ 1 ~B' u y u , which indicates equivalency to the estimator f yk 
obtained with penalized B-spline fitting given by (22) for pL' v ALu = AD g+ i (see Claeskens 
et al. (2009) for the expression of Lu satisfying this equation). 

In a design-based approach, Claeskens et al. (2005) proved that the NMA estimator t y> BS is 
the population total of the design-based predictions f y ^ = c '( z k)(C' s II~ 1 C s +pA)~ 1 C' s II~ 1 ~y s . 
They also proved that t yiBS fulfils properties (12) and (13). 



4.2 Asymptotic properties of the B-spline estimator of totals un- 
der the sampling design 

In the following, we study the asymptotic properties of t Vt Bs under the sampling design. We 
first provide a lemma concerning the convergence of y ,\- The proofs are based on the results 
provided by Goga (2005) for the unpenalized B-spline estimator and on the fact that the 
inverse of the matrix ^B^Bj/ + is of order O(K) for the penalized B-spline estimator 
(lemma 1 from Claeskens et al, 2009). 

Lemma 5. (a) Assume assumptions (AJ^)-(b) and (Bl), (B2)-(a) and (B3) from the Ap- 
pendix. Then, H^aII = O^K 1 / 2 ). 

(b) Assume assumptions (A3), (AJ^)-(b), (A5) and (B1)-(B3) from the Appendix. Then, 

E P (\\o y , x -~e yA \\ 2 )=o(^y 

where \ \ ■ \ \ is the usual euclidian norm. 

The proof is provided in the Appendix. We note that for B-spline functions of order m = 1 
and A = 0, we obtain a poststratified estimator with a number of poststrata going to infinity. 
In this context, lemma |5j (b) provides a detailed theoretical justification for the poststrati- 
fication example in Deville (1999, p. 196). We note also that to obtain the convergence of 
iy,BS> Claeskens et al. (2005) assume that the result from lemma [5J (b) holds. Finally, we 
note that GREG estimators may be viewed as a special case when the number of knots is 
fixed. Papers dealing with this issue usually assume that the regression coefficient satisfies 
the results from the above lemma (see for example Robinson and Sarndal, 1983, or Isaki and 
Fuller, 1983). A similar result was proved by Cardot et al. (2012). 
Using lemma [5j we derive the following results. 

Proposition 6. Assume assumptions (A3), (AJ^)-(b), (A5) and (B1)-(B3) from the Ap- 
pendix. Then, 
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(a) E p \jt(i y , BS -t y )\=0{{K/nyi*)). 

(b) i (i y ,BS -ty) = ^ (W - *») + °> _1/2 ) ufcere = £ s ^ + £t/ /»,*• 

The proof is provided in the Appendix. Using the Markov inequality, we see from the first 
point of proposition [6] that i Vi Bs is asymptotically design-unbiased for t y and i/ra-consistent 
as (i y ,BS — t y )/N = O p (n~ l l 2 ). The second point provides an approximation of i y> BS by the 
nonparametric generalized difference estimator t* diff . 



4.3 Calibration with penalized splines 



The spline approach has some interesting calibration properties. Under the unpenalized 



B-spline framework, the weights w ks given by (26) satisfy the calibration equation for the 



known population total of B-spline functions, namely 

^2 w ks Bj(z k ) = ^ B j( z k), for all j = 1, . . . , q. 



u 



This relation is easily obtained using (20) (Goga, 2005). Because the spline space S^m is 
spanned by the B-spline functions Bj, these weights will be calibrated to the total of any 
polynomial z r of degree r < q = K + m. In particular, E s w ks = N and E s w ks z k = Ylu Zk - 
Claeskens et al. (2005) prove that using the penalized splines and the truncated polynomial 
basis functions 1 provides weights that are also calibrated for the finite population totals of 
the polynomial basis functions 1, z, z 2 , . . . , z m ~ x . 



4.4 Nonparametric penalized spline estimation for nonlinear pa- 
rameters 

We now consider the nonlinear parameter $ estimated by §bs — T(Mbs) with Mbs — 
E s w ksby k and the weights w ks given by ( 
two-step procedure. The first-step linearization is given in theorem [3] provided that the 
assumptions (A*) and (A**) from lemma ^ are fulfilled. These assumptions are crucial 
because they ensure the convergence of some nonparametric estimator of M to the true mea- 
sure M according to the distance dt v . Using classical assumptions from a B-spline framework 
(assumptions (B1)-(B3) from the Appendix) and mild assumptions regarding the sampling 
design (assumptions (A3) and (A5) from the Appendix), we prove in theorem [7] below that 
(A*) and (A**) are verified. The proof is basically based on lemma [5] and the fact that the 
distance d tv is defined for uniformly bounded functions h 6 H, ensuring that the assumption 
(A4)-(b) is automatically fulfilled. 



25). As in section s to linearize $bs we use 



18 



By conducting this first linearization step, we see that the nonparametric B-spline esti- 
mator <& B s will be approximated by the nonparametric B-spline estimator of the total of the 
linearized variables Uk given by 

,* / Uk ~ 9 u ,k . * 

t u ,Bs = w s u s = > + > g Uik , 

s nk U 

where g* u>k = b'(z k )d u , x with U , X = (B'Jl^B, + AD p ) -1 B' s Il7 1 u s . 

The second-step linearization consists of providing an approximation of t* BS by a non- 
parametric generalized difference estimator, 

u k — g u ,k 



t u,dis — 7, r > g u ,k 



^k 

o u 

where g u , k = h'{zk){B' v Bu + AD p )~ 1 B' [7 uu- To obtain this result, we state in theorem [7j 
(b) a supplementary assumption regarding the linearized variable u k . Goga and Ruiz-Gazen 
(2012) prove that the linearized variable u k of the odds-ratio satisfies this assumption. 

Theorem 7. Suppose that the sampling design satisfies assumptions (A3) and (A5). In 
addition, assume that (B1)-(B3) hold. 

(a) Assumptions (A*) and (A**) from lemma^ are fulfilled. 

As a consequence, d tv (M B s/N, M/N) = O v (n~ l l 2 ). Moreover, if the functional T sat- 
isfies (Al) and (A2), then N- a ($ BS - $) = N- a (i U)BS - t u ) + o p {n-^ 2 ). 

(b) Suppose that the linearized variables are such that for all k G U, N~ a+1 u k satisfy 
(A4)-(b). Then, N~ a (t u , BS - t u ) = N~ a (i u , m -t u ) + o^n- 1 ' 2 ). 

The proof is provided in the Appendix. 

5 Variance estimation 

In this section we undertake a detailed study of the variance estimation of $ np . We first give 
the functional form of the variance of t Vt HT a s well as of its variance estimator and we propose 
a variance estimator for <3> np and assumptions under which this estimator is consistent. 

The Horvitz-Thompson variance V p (i yyHT ) = Y^uY^u ^ki{Vk/^k){yi/n) for A H = tt h - 
7Tfc7T; is a quadratic form that can be written as a functional of some finite and discrete 
measure. We can write the variance as follows (Liu and Thompson, 1983), 

V p {ty iHT )= J2 ^(Vkiyi) (27) 
(k,i)eu* 
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where U* = {(k,l), k,l = 1, . . . N} and i/;(y k ,yi) = A k i(y k /n k )(yi/ni) is a bilinear function of 
y k and y\. It follows from (27), that the Horvitz-Thompson variance V p is the finite population 



total of ip(y k ,yi) over the derived synthetic population U* of size N* = A^ 2 . This variance 
can be put in a functional form as follows 

V p (i y>HT ) = T*(M*) = J iP{y,y)dM*(y,y) 

where M* = 5^(jfej)ei7* ^(yk,yi)- Note that T* is a functional of degree 1 with respect to M*, 
namely T*(M*/N*) = T*(M*)/N*. A sample in this population U* is s* = {(k,l),k,l G s} 
and has size n* = n 2 . Moreover, the first-order inclusion probabilities over the synthetic 
population U* are 7r* fc ^ = ir k i, which are exactly the second-order inclusion probabilities 
with respect to the initial sampling design p(s). The measure M* is estimated on s* by 
M * = £(*,!)£*• s (vh,vi)/*u = £ s * w \ki) 5 {y k m) where w *m) = 1 / 7r ki- The resulting estimator of 
Vp(t y ,HT) is as follows 



V p (i y ,HT) = T*(M*) = j ^(y,y)dM*(y,y) = £ 



Afc/ y k yi 



(fe,t)6s* 

This is exactly the Horvitz-Thompson variance estimator, as J2( k z)e s * * s ec L ua l to J2kes J2ies ■ 
Moreover, the functional T* is Frechet differentiable, with first derivative given by IT*(M*, y) ■ 

1>(v,y)- 

Consider now the asymptotic variance AV p {Q np ) of $ np given by 

AV p @ np ) ^VVA ^^ (28) 

z — ' z — ' 7T k 7T; 

U U 

where u k is the linearized variable of $ and g Ujk = q,' k uu for uu = (u k ) k eu- We recognize 
the Horvitz-Thompson variance of the total of the population residuals e ks = u k — g Ujk . We 
suggest estimating the variance of $ np by using the Horvitz-Thompson variance estimator 
with u k replaced by the sample estimators u k , 

hk * ~ ^ (29) 

^ ^ TC ki TC k Hi 

s s 

where cju,k = qL^s * s the sample estimate of g Ujk = q. k uu- The Horvitz-Thompson variance 
estimator with true linearized variables given by 

av p (K) = yy^i Uk ~ Ul ~ hi . ( 3o) 

^ ~ T^kl TTfc T\\ 
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The three expressions of variance above depend on the population fits residuals ek S} for all 
k 6 U. It follows that we may write AV p (<& np ) as a functional of M* depending on parameter 
&u = (eks)keu, 

AV p ($ np ) = T*(M*,eu). 

Furthermore, the Horvitz-Thompson estimator AV p (§ np ) and the variance estimator V p (Q np ) 
can be treated in a functional form as follows 

AV p ($ np ) = T*(M*, eu ), V p ($ np ) = T*(M*,e v ). 

Note that e v = (e fcs ) fcgC / is the vector of sample-based fit residuals with e^s = Uk — 9u,k, 
for all k G U. Theorem 3 from Goga et al. (2009) allows us to establish under additional 



assumptions that the variance estimator (29) is n-consistent for the asymptotic variance 



Theorem 8. Assume that assumptions (A3) and (A5) from the Appendix hold. Also as- 
sume that N 1 ~ a eks = 0(1) holds uniformly in k and nN~ 2aS ^^(ek s — e-ks) 2 — o p (l). If the 



v 



Horvitz-Thompson variance estimator AV p (<& np ) is n- consistent for AV p (Q np ), then the vari- 
ance estimator V p {^ np ) is also n-consistent for AV p {Q np ) in the sense that nN~ 2a {V p {^ np ) — 
AV p ($ np )) = o p (l). 

The proof is given in the Appendix. Note that because the functional T* is Frechet differ- 
entiable, the n-consistency of the Horvitz-Thompson estimator AV p ($> np ) for AV p (§ np ) may 
also be derived with assumptions on fourth moment of and on fourth-order inclusion 
probabilities. The reader is referred to Breidt and Opsomer (2000) for additional details. 

6 Empirical results 

Let us consider a data set from the French Labor Force surveys of 1999 and 2000 as the 
finite populations of interest. The data consist of the monthly wages (in euros) of 19,378 
wage-earners who were sampled in both years. The study variable (resp. the auxiliary 
variable Xk) is the wage of person k in 2000 (resp. 1999). The objective of the simulation 
studies is to investigate the finite-sample performance of the regression spline estimators for 
two nonlinear parameters of interest and two different survey designs. We concentrate in 
practice on the simple approach of regression B-splines and do not consider the penalized B- 
splines with A > 0. The empirical study of penalized splines raises the problem of estimating 
the parameter A which is beyond the scope of the present paper. We illustrate the efficiency 
of the regression B-splines estimators compared to other estimators, and we also confirm the 
possibility of conducting valid inference using variance estimators as detailed in the previous 
section. 
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The parameters to estimate include the mean, the Gini index and the poverty rate for 
the wages in 2000 using the wages in 1999 as auxiliary information. The poverty rate is the 
proportion of individuals whose wages are below the threshold of 60% of the median wage 
and correspond to the low-income proportion studied in Berger and Skinner (2003). The Gini 
index and the low-income proportion are the complex parameters to be estimated and we 
provide results for the mean as a benchmark. Note that details on the low-income proportion 
estimator and its associated linearized variable can be found in Berger and Skinner (2003) 
and are not provided in the present paper. In subsection 6T we focus on simple random 
sampling without replacement and in subsection 6.2, we focus on a stratified simple random 
sampling without replacement. We consider the following estimators for each parameter: 

- the Horvitz-Thompson estimator (HT), which does not incorporate any auxiliary informa- 
tion, 

- poststratified estimators (POST) with a different number of strata bounded at the empir- 
ical quantiles for 1999 wages, 

- the GREG estimator (GREG), which takes into account the 1999 wages as auxiliary infor- 
mation using a simple linear model, 

- B-spline estimators (BS(m) where m denotes the spline order), which take into account 
the wages from 1999 as auxiliary information by using a nonparametric model with different 
numbers of knots (K) located at the quantiles of the empirical distribution for wages from 
1999. The m = 2 and m = 3 orders are considered. 

The poststratified estimator is an example of a B-spline estimator with order m = 1. 
The number of strata correspond to the number of interior knots K plus one. 

To use the regression B-spline estimators we propose in a complex survey, and derive 
confidence intervals, the user must be able to calculate the weights given in equation (26) 
and the residuals itk — gu,k of equation (29 ). The weights depend on a spline basis that is easy 
to obtain using for instance the transreg procedure in the SAS software or the functions 
spline. des or bs from the splines package in the R software. Then, it is possible to use 
standard calibration algorithms by simply providing the (m + K) B-spline basis functions as 
auxiliary variables for calculating the calibrated weights that correspond to equation (26). 
These weights are needed to calculate the substitution estimator of the parameter of interest 
(e.g. the expression (15) for the Gini index). To estimate the variance, the linearized vari- 
ables associated with the parameter have to be estimated. For several inequality indicators, 
including the Gini index and the low-income proportion, some SAS macro programs are 
freely available on the web site of Xavier d'Haultfceuille. Similar functions are available in 
the R language upon request from the authors of the present paper. Once the linearized vari- 
able is estimated, the residuals of this variable against the auxiliary variable using regression 
splines are calculated; this can be accomplished with the transreg procedure in the SAS 
software. Then, by using the residuals as if they were the study variable in standard vari- 
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ance estimation tools for complex surveys, the user can obtain the estimated approximative 
variance and derive confidence intervals. 

For each simulation scheme, we draw NS samples according to the sampling design 
and compare the finite-sample properties of the HT estimator, the GREG estimator, the 
POST and the BS(2) and BS(3) estimators. We set knots at the quantiles of the empirical 
distribution of the auxiliary variable in the sample. We also compared the results with knots 
set at the quantiles of the empirical distribution of the auxiliary variable over the entire 
population. Both results are very similar; thus, we report only on the first method. For the 
POST, BS(2) and BS(3) estimators we tried different numbers of knots K but only report 
the results for K = 2 and K = 4. Note that in the tables, the results for K = 2 and K = 4 
are reported in the same columns and separated by a dash. For the poststratified estimator, 
K = 2 (resp. K = 4) corresponds to 3 (resp. 5) strata. To summarize, in the following, 
we compare eight estimators (HT, GREG and POST, BS(2) and BS(3) with K = 2 and 
K = A). 

There are several ways to estimate the linearized variable (see section [5J. In this section, 
the results are almost the same, regardless of whether we use the simple HT weights, the 
GREG weights or the B-spline weights for estimating the linearized variable. We recommend 
using the simplest weights (that is, the HT weights), which is what we do in the present 
study. 

Estimators performance of 9 for a parameter 9 is evaluated using the following Monte- 
Carlo measures: 

100 NS 

• Relative bias in percentage: RB = x ^^(#i — 9)/9. 

' i=l 

• Ratio of root mean squared errors in percentage: 



RRMSE = 100 x 



NS 



NS 



\ 1=1 \ 1=1 

Monte-Carlo Coverage probabilities for a nominal coverage probability of 95%. 



6.1 Simple random sampling without replacement 

The first survey design we consider is simple random sampling without replacement with 
three sample sizes (n = 200, n = 500 and n = 1000). The number of simulations is 
NS =3,000. The eight estimators are compared and relative biases and ratios of the roots 
of the mean squared errors are provided in Table [T] for the different parameters and sample 
sizes. 
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Not surprisingly, for complex parameters, the largest efficiency gain is observed when 
the B-spline estimators are compared to the HT estimator without auxiliary information. 
Because the wages from 2000 are almost linearly related to the wages from 1999, considering 
the B-spline estimator instead of the GREG estimator does not improve the performance of 
the mean estimation. However, regarding the Gini index and the low-income proportion, the 
incorporation of auxiliary information using GREG estimators does not improve efficiency 
compared to the HT estimator while using a B-spline approach improves the results especially 
for spline functions of order m — 2. When comparing the POST estimator with the BS(2) 
and BS(3) estimators we notice that there is quite a large gain in efficiency when order m = 2 
is used instead of m = 1, while there is an efficiency loss when m = 3 is used instead of 
m = 2, especially for sample sizes smaller than 1,000. Moreover, for m = 2 and m = 3, the 
results do not depend heavily on the number of knots and are similar for K between 2 and 
4 while for the poststratified estimator, there are large variations in the results, regardless 
of whether we consider 3 or 5 strata. The coverage probabilities in table [2] illustrate that 
valid inference can be carried out using B-spline estimators as long as the spline order is 
not too high, especially when the sample size is not very large. No problems are detected 
for B-splines of order m = 1 and order m = 2 even when the sample size is n = 200; 
however for m = 3 and n = 200, the coverage probabilities for the Gini index estimation 
are approximately 75% which is quite far from the 95% nominal probability. This result 
indicates that for a moderate sample size, the variance may be underestimated when the 
order of the splines is larger than two. The results are not given for m = 4 but we have 
observed that the problem worsens when we increase the order of the splines. This is not 
really surprising due to double linearization and nonparametric estimation. 

6.2 Stratified simple random sampling without replacement 

For each simulation, we draw NS =5,000 samples from the French Labor Force population 
according to a stratified simple random sampling design without replacement. We compare 
the finite-sample properties of the eight estimators considered in the previous subsection. 
The strata are spatial divisions of the French territory in six "regions" that correspond to 
the major socio-economic regions of metropolitan France as defined by Eurostat. These 
regions are the first level of the nomenclature of territorial units for statistics classification 
(NUTS-1). For our example, we grouped the Northern and Eastern regions together and we 
grouped the Mediterranean and the Southwestern regions together. The sample size inside 
each stratum is 200 making the total sample size 1200. Thus, we used an unequal probability 
design with a sample rate inside the strata that varied from 5 to 9.3%. 

As previously described, we set the knots at the quantiles of the empirical distribution of 
the auxiliary variable in the sample and we estimate the linearized variables using the HT 
weights. The simulation results are reported in Table [3] and [4] and the conclusions are similar 
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Table 1: RRMSE (RB) of HT, GREG and POST, BS(2) and BS(3) with K = 2 - 4 for the 
mean, the Gini index and the low-income proportion 



Parameter 


n 


HT 


GREG 


POST 


BS(2) 




BS(3) 


Mean 


zUU 


1UU ) 


OO ) 


71 (r\\ (r\\ 


38 (0) 


- 37 (0) 




(C\\ /II (r\\ 

[(J) - 41 [0) 




500 


100 (0) 


40 (0) 


73 (0) - 65 (0) 


40 (0) 


- 39 (0) 


38 


(0) - 39 (0) 




1,000 


100 (0) 


40 (0) 


73 (0) - 66 (0) 


40 (0) 


- 40 (0) 


38 


(0) - 39 (0) 


Gini index 


200 


100 (1) 


96 (1) 


92 (1) - 80 (1) 


53 (2) 


- 53 (2) 


70 


(3) - 70 (3) 




500 


100 (1) 


93 (0) 


93 (1) - 85 (1) 


50 (1) 


- 50 (1) 


59 


(1) - 56 (1) 




1,000 


100 (0) 


92 (0) 


93 (0) - 86 (0) 


49 (0) 


- 48 (0) 


55 


(1)-51 (1) 


Poverty rate 


200 


100 (2) 


95 (0) 


92 (0) - 80 (0) 


65 (1) 


- 65 (1) 


72 


(1) - 63 (1) 




500 


100 (0) 


95 (0) 


88 (0) - 78 (0) 


64 (0) 


- 64 (0) 


68 


(0) - 62 (0) 




1,000 


100 (1) 


94 (0) 


89 (0) - 78 (0) 


64 (0) 


- 64 (0) 


67 


(0) - 61 (0) 



Table 2: Coverage probabilities (in %) for HT, GREG and POST, BS(2) and BS(3) with 
K = 2 - 4 



Parameter 


n 


HT 


GREG 


POST 


BSPL(2) 


BSPL(3) 


Mean 


200 


94 


95 


93 - 


92 


93 - 


93 


90 - 


88 




500 


95 


94 


93 - 


94 


93 - 


93 


91 - 


91 




1,000 


95 


95 


94 - 


93 


94 - 


94 


93 - 


93 


Gini index 


200 


94 


93 


94 - 


94 


89 - 


87 


74 - 


75 




500 


93 


93 


93 - 


94 


91 - 


90 


83 - 


85 




1,000 


95 


94 


95 - 


94 


94 - 


93 


88 - 


90 


Poverty rate 


200 


94 


95 


95 - 


95 


95 - 


94 


94 - 


94 




500 


93 


95 


95 - 


94 


95 - 


95 


96 - 


95 




1,000 


94 


95 


96 - 


96 


95 - 


96 


96 - 


95 
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Table 3: RRMSE (RB) of HT, GREG and POST, BS(2) and BS(3) with K = 2 - 4 



Parameter 


HT 


GREG 


POST 


BS(2) 


BS(3) 


Mean 


100 (0) 


40 (0) 


73 (0) - 66 (0) 


40 (0) - 40 (0) 


40 (0) - 40 (0) 


Gini index 


100 (0) 


93 (0) 


94 (0) - 88 (0) 


50 (0) - 50 (0) 


55 (1) - 52 (1) 


Poverty rate 


100 (0) 


93 (0) 


88 (0) - 77 (0) 


65 (0) - 64 (0) 


68 (0) - 62 (0) 



Table 4: Coverage probabilities for HT, GREG and POST, BS(2) and BS(3) with K = 2 - 

4 



Parameter 


HT 


GREG 


POST 


BS(2) 


BS(3) 


Mean 


95 


95 


95 - 94 


94 - 94 


93 - 92 


Gini index 


95 


95 


95 - 95 


93 - 93 


89 - 90 


Poverty rate 


94 


95 


95 - 95 


95 - 95 


96 - 95 



to those obtained from the simple random sampling design without replacement when the 
size of the sample is n = 200 which corresponds to the sample sizes inside each stratum. 
It is beneficial to use the available auxiliary information when estimating the mean but 
there is no need to use nonparametric estimators because they are not more efficient than 
the GREG estimator. However, for complex parameters, using a GREG estimator to take 
auxiliary information into account is not worthwhile in terms of variance while important 
gains can be made by using B-spline estimators. The empirical coverage probabilities are all 
very good except for the Gini index B-spline estimator of order three with values equal to 
89-90% which confirms the problem of variance underestimation for moderate sample sizes 
and splines of order three. 

Based on this example we do not recommend using high order values for B-spline regression, 
especially when the sample sizes are smaller than 500. However, choosing m = 2 instead of 
m = 1 (which corresponds to poststratification) leads to a clear improvement in terms of 
efficiency for complex parameters such as the Gini index or the low-income proportion, and 
we recommend this choice. 



7 Discussion 

In this paper we considered the important problem of nonlinear parameter estimation in 
a finite population framework by taking into account the survey design and a unique aux- 
iliary variable known for all the population units. Examples of nonlinear parameters are 
concentration and inequality measures, such as the Gini index or the low-income proportion. 
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We proposed a general class of substitution estimators that allows us to take into account 
the auxiliary information via a nonparametric model-assisted approach. The asymptotic 
variance of this class of estimators was derived, based on broad assumptions, and variance 
estimators were proposed. Our main result was that the asymptotic variance depends on the 
extent to which the auxiliary variable Zk explains the variation in the linearized variable Uk- 
Because linearized variables of nonlinear parameter are likely to be nonlinearly related to 
auxiliary information, a nonparametric approach is recommended. The proposed estimators 
are based on weights that are flexible enough to increase the efficiency of finite population 
totals estimators for any study variable and to allow the consideration of parameters that 
are more complex than totals. Moreover, the penalized B-spline estimators were studied in 
detail, and the theoretical results were confirmed for regression B-spline estimators using 
one case study. 

Our proposal can be extended in several different ways. In particular, further research 
can extend this proposal to include multivariate auxiliary information by means of additive 
models, as in Breidt et al. (2005), or single index models as in Wang (2009). 

Acknowledgement: we are grateful to Patrick Gabriel for his precious help for lemma 
1, to Herve Cardot for helpful discussions and to Didier Gazen for his assistance with the 
simulations. 

Appendix: assumptions and proofs 

Assumptions on functional T and on sampling design. 

(Al). The functional T is homogeneous, in that there exists a real number a > 0, dependent on T 
such that T(rM) = r a T(M) for any real r > 0. We assume also that limjv-».oo N~ a T(M) < 
oo. 

(A2). The functional T is Frechet differentiate at M/N; that is, there exists a functional T(M/N; A) 



that is linear in A such that 
withd(£,#) — >0. 



G\ (M\ „ M G-M 
T 



N \N U N 



o d 



G M 

AT' iV 



We note that the strong assumption of Frechet differentiability can be weakened to compact or 
Hadamard differentiability. However, for Hadamard differentiability, functionals are considered 
with respect to the empirical distribution function and the distance a\ Y should be replaced by 
the sup norm. Supplementary assumptions need to be supposed in order to have the consistency 
of the estimator of the empirical distribution function. Motoyama and Takahashi (2008) study 
the asymptotic behavior of Hadamard statistical functionals but only for simple random sampling 
without replacement. 



27 



(A3) 



vre (0,1). 



(A4). (a) linx/V 1 Y^uVk < 00 w ^ n ^-probability 1. 

(b) sup fe6 [/Z/fc < C with C a positive constant not depending on N. 

(A5). min7Tfc > A, minj keu ^ik > A* with A, A* with some positive constants and limjv_ 5 . 00 n max Wik— 

k£U ' " i^=k£U 

TTiir k \ < 00. 

Assumption (A3) and (A5) deal with first and second order inclusion probabilities and are rather 
classical in survey sampling theory (see also Robinson and Sarndal, 1983 and Breidt and Opsomer, 
2000). They are satisfied for many sampling designs. Assumption (A4)-(a) is a regularity condition 
necessary to get the consistency results. Some results need the stronger assumption (A4)-(b). 

Assumptions on B-splines 

(Bl). There exists a distribution function Q(z) with strictly positive density on [0,1] such that 
sup zg [ 0j i] \Qn{z) — Q{z)\ = o(K~ l ), with Qn(z) the empirical distribution of [zi)f =1 . 



(B2). (a) K = o(N); 

(b) K = 0(n a ) with a < 1/3. 

(B3). K p = (K + m-p)(Xc) 1 /^N- 1 / < - 2 P^ < 1 where c = c(l + o(l)) with c a constant that depends 



only on p and the design density. 

These assumptions are classical in nonpar ametric regression (Agarwal and Studden, 1980, Bur- 
man, 1991, Zhou et al, 1998); (Bl) means that asymptotically, there is no sub-interval in [0,1] 
without points Z{ and (B2) ensures that the dimension of the B-spline basis goes to infinity but not 
too fast when the population and the sample sizes go to infinity. Assumption (B3) concerns the 
penalty A as used by Claeskens et al. (2009). 

Proofs of results from section [3] 

Proof of Lemma [TJ Now, let h S H and let = l{fc £s } be the sample membership. Following 
the same lines as in Breidt and Opsomer (2000), we have, 




hdM/N 



2 




< 



1 - A wmax fc ^|A M 



XN X 2 n 
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uniformly in h by assumption (A3),(A5) and using the fact that h £ H. 



Proof of Lemma [2l We have 




From the proof of lemma [IJ we see that the first term from the right-side is of order O^" 1 ^ 2 ) 
uniformly in h because (1/N) Y^u( x k ~ fx,k) 2 < (2/AT) Y^u( x l + fx,k) — 2 (1 + C) by construction 
of Xk and assumption (A*). The result follows because | J hdM np /N — J hdM/N\ = O p (n -1 / 2 ) 
uniformly in h G 



Proof of Theorem [3] Under assumption (A2), we provide a first-order von-Mises (1947) expansion 
of T in M np /N around M /JV, 

Using the fact that for a functional of degree a (assumption Al), we have IT(^,£) = iV 1-a • 
JT(M,f) (Deville, 1999), we write 

N- a T{M np ) = N- a T{M) + N- a J IT (M, £) d(M np — M)(£) + o p (ra -1 / 2 ) (31) 



since yM np /N,M/Nj = O p (n x / 2 ). Now, = IT(M,yk) and hence, relation (31) becomes, 

Proofs of results from section |4] 

We state below several lemmas useful for the proofs of our main results. For a matrix A 

( a ij)ij=ii we consider the norm defined by ||A||oo = max^ =1 Y^j=i \ a ij\ an d the trace norm ||A||2 
(trace(A'A)) 1 /2. 

We denote by Ha = ^B'^B^ + ^D p and by Ha = ^B' s nj 1 B s + ^D p its estimator. 

Lemma 9. Assume assumptions (Bl), (B2)-(a) and (B3). Then, 

1. 1 1 ^(B'^Bf/)! |oo = 0(K~ l ), (lemma 6.3 from Agarwal and Studden, 1980). 
We also have \ \{j f & u B u )- 1 \\ 00 = 0{K), (lemma 6.3 from Zhou et al, 1998). 
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2. IK^B^Bp + ^Dp)" 1 !^ = HH^Iloo = 0(K) (lemma 1 from Claeskens et al, 2009) 
Lemma 10. (Goga, 2005) Assume (A3), (A4)-(a), (A5) and (Bl), (B2)-(a). Then, 

1. J B p ||^(B' s n; 1 B s -B' [/ B l/ )||| = 0(i) 

2. J B p ||^(B' s n- 1 y s -B' f/ y c/ )|| 2 = 0(i) 

Proof of lemma [5] When yf. is uniformly bounded (assumption A4,b), we have, using lemma 3 
(a) (Goga, 2005) that 

ii^ y Hi 2 <^ii£b ( ^ii 2 <^ (32) 

u 

since for k,l G U with |A; — l\ > m we have Bj(xk)Bj(xi) = 0. 
For (a), y> \ is bounded following Goga (2005), 

fell < HH^iU-iKi/^B^H 

= OiK- 1 ' 2 ) (33) 



by lemma |9j-(b) and relation (32). Furthermore, we have 
\\&y,\ — 6y,\\\ 2 

< IIH- 1 - H- 1 ^ • n^B^njV.II 2 + Willie -11^ (BinrV. - b^) || 2 (34) 

Under the assumption (A4)-(b), H^B'^njVsll 2 is bounded by \\jjJ2u h ( z k)\\ 2 = 0(K- V ). We 
have that 



-H- 1 (1 (B'^^B, - B'uBu)) (h- 1 ! (B'^^B, - B^) + I,) H 



-l 

-l 

A 



and HH" 1 ^ (B^n;:^ - B^Bu) = o p (l) for K = 0{n a ) with a < 1/3, implying that 
|| (H^ 1 ^ (B^n^^s - B'vBu) + Iq)" 1 | |oo < 1. Using lemma 10 -(a), we obtain that 



From lemmas [9] and 10 we obtain that 

e p (WhL • ||1 (B' s n r V, - B' uyu ) || 2 ) = o (^) 
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Finally, we have that 



Ep\\O y ,\ - 0y,\\\ -O 



n J 



Proof of proposition [6j Consider first (b). Using the same lines as in the proof of lemma [I] and 
the fact that ||b(zfe)|| < 1 for all k G U (Burman, 1991), we obtain that 



0(n 



-l/2\ 



(35) 



Furthermore, 
1 



N 



ty,BS — ty,diff) 



< E n 



• \ \0y,\ — 0y,\\ 



< 



En 



1 h'(z k ) 



" Ep\\9y : \ — 0y,\\ 



o 



^3/2 

V n 

oin- 1 ' 2 ) 



(1=\ . 0((K 3 /n) 1/2 ) = (-^=) ■ 0((n 3a - 1 ) 1 / 2 ) 



by (35) and lemma [5J- (b). Then, the result follows by using the Markov inequality, 
(a) Now, we consider the error t y bs ~ ty We write 



N 



{i y ,BS ~ t y ) 



N 



{iy,HT ~ t y ) 



jf E b '(^) (| - l ) - »».a) - ^ E b '(^) ( 



4 



1 



By assumptions (A3), (A4-b) and (A5), we have that E p \j^ [iy,HT — ty) I = 0(n 1 ^ 2 ). Moreover, 
using relation (35) and lemma jsj (a) we have E p |^ Y.u {k ~ l) b'^fe)^! = O^/n) 1 / 2 ) which 
implies that 



En 



± (iy,BS -t y )\< 0{n^ 2 ) + OiK^n- 1 ) + O^n" 1 / 2 ) = 0({K/nfl 2 ) 



by the fact that (K/n) 1 ' 2 > n l / 2 > K 3 / 2 n 1 using assumption (B2). 

Proof of Theorem [7j (a) We check that assumptions (A*) and (^4**) are fulfilled. We have 
6 XiX = H^ 1 (X; c/ b(z fc )x fc /iV) with \x k \ = \h{y k )\ < 1 for all k € U. Following g and Q, we 
obtain that ||0 x ,a|| = 0(iT 1//2 ) uniformly in /i and 



Jj^2fx,k = Jrdx,X S 'u' B udx,>> < ll^,Al| 2 ||^y B C7 B t/||oo = 0(1 

u 



(36) 
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uniformly in h. Now, we check the assumption {A**), namely E p ^ ( ^ — lj (f x ,k — f x ,k] 
0(n -1 / 2 ) uniformly in h. 
We have 



En 



1 (fx,k — fx,k) 



< E n 



N ^ 



h 



1 b'fo) 



< 



1 



a, 



l -y 



h 



N \ir k 



i b'(z fc ; 



The first term from the right-side does not depend on h and is of order 0(n _1 ) (equation (35 )). For 
the second term from the right-side, we can use the proof of lemma ([5]), more exactly the equation 
(34), and the fact that sup fcg; y < 1 to obtain 



E P \\0 X \ — X 



'K 3 \ 

O ( j uniformly in h. 



Finally, we obtain that E p 



hT.u(k- l ) (/*.* " /»>*) = o(n- 1 / 2 ) for K = 0(n a ) with a < 1/3. 



(b) We write equation ([7]) as follows: 

N- a (t* u , BS ~ t;, diff ) = N- a Y, (~ - - Kx) = o P (n- 1/2 ) 



because N 1 ^2 U 
by lemma [5j 



4 



1^ h'{z k ) = O p (n- 1 / 2 ) (equation (35)) and N- a+1 {0 u -~0 u ) = O p {K^ 2 rT l l 2 ) 



Proof of Theorem [8| The proof follows the same basic steps as in Theorem 3 from Goga et al. 
(2009) and result 4 from Chaouch and Goga (2010). Let 



A N = Vp($ np ) - AV p ($ np ), B N = AV p ($ np ) - AV P {<S> 



up j 



written as 



with AV p {<& np ) given by (30) and let also c k i = — ^._^L_i_ Furthermore, the quantity A^ can be 



A 



N 



^ ^ Cki(ek s eis - e ks eis) 
u u 

^ ^ Cki(eks - e ks )(ei s - e is ) + 2 ^ ^ c fc ;(e fcs - efc s )e is 



u u 

A±N + A 2 N 



u u 



Now, 



n 



Ai N \ < 



1-A n ^ 2 nmax|A H |^ 2 



A 2 A^ 2c 



u 
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by assumptions (A3) and (A5). Using the same arguments as above, we obtain nN 2ol \A2n\ = 
o p (n _1 ). Hence, nN~ 2a \A^\ = o p {n~ l ) and the result then follows because uN~ 2o! Bn = o p (l) 



N 2 



-(V p ($ np ) - AV p ($ np )) < j^(\A N \ + \B N 



n 
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